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ABSTRACT 

Using highly- accurate A^-body simulations, we explore the evolution of multiple massive black holes 
(hereafter, MBHs) in a primordial galaxy that is composed of stars and MBHs. The evolution is 
pursed with a fourth-order Hermitc scheme, where not only three-body interaction of MBHs but also 
dynamical friction by stars are incorporated. Initially, ten MBHs with equal mass of 10'' Mq are set 
in a host galaxy with 10^^ Mq. It is found that 4-6 MBHs merge successively within 1 Gyr, emitting 
gravitational wave radiation. The key process for the successive merger of MBHs is the dynamical 
friction by field stars, which enhances three-body interactions of MBHs when they enter the central 
regions of the galaxy. The heaviest MBH always composes a close binary at the galactic center, which 
shrinks owing to the angular momentum transfer by the third MBH and eventually merges. The 
angular momentum transfer by the third MBH is due to the sling-shot mechanism. We find that the 
secular Kozai mechanism does not work for a binary to merge if we include the relativistic pericenter 
shift. The simulations show that a multiple MBH system can produce a heavier MBH at the galactic 
center purely through A^-body process. This merger path can be of great significance for the growth of 
MBHs in a primordial galaxy. The merger of multiple MBHs may be a potential source of gravitational 
waves for the Laser Interferometer Space Antenna (LISA) and pulsar timing. 

Subject headings: black hole physics — galaxies: nuclei — methods; numerical 



1. INTRODUCTION 

In the central regions of galaxies, massive black holes 
with 10^ - IO^Mq (hereafter, MBHs) are found to reside. 
They are believed to coevolve with their host galaxies, 
since the masses of MBHs correlate with the mass and ve- 
locity dis persion of the sp heroidal components of the host 
galax i es (iKormendv fc R ichstonc 1995; Mago rrTan et al.l 
19981; iFerrarese fc Merritt .2000; .Tremaine et al.l 120021; 
Marconi fc HuntI l2003f l Such MBHs seem to acquire 



most of their masses through gas accre tion process in 
the fi nal evolutionary stage (Soltan 1982; Yu fc Tremaind 
2002'). But, the growth from an order of magnitude 
smaller black holes is still unresolved. 

In the last decade, man y quasars at re dshift z ~ 6 
have been observed (e.g.. Iran et"all[200l . This sug- 
gests that MBHs with 10^ Mq have formed when the age 
of the Universe is only 1 Gyr. If seed blac k holes are 
the remnants of fi rst stars with ~ lOOMfn (lAbe et al 
2000; Nakamura fc_ Umemural 120011 ; iBromm et al.l 1200 



i 



Yoshida et al the seed black holes cannot grow to 



10 M0 with continuous Eddington accretion rate over 1 
Gyr. One solu tion is the growth by super-Eddington 
accre tion (e.g., lAbramowicz et all 119881 ; lOhsuga et al.l 
|2005[ ). However, gas accretion onto the seeds should be 
intermittent, and o n average could be lower th an the Ed- 
dington accretion (jMilosavlievic et aLll2009al lbl). 

Hence, it is worth considering the possibility of black 
hole mergers for the growth of MBHs. In the cold dark 
matter cosmology, a massive galaxy forms through the 
multiple merger of subgalaxies. If subgalaxies possess 
MBHs, a massive galaxy should contain multiple MBHs 
shortly after the merger. On the ot her hand, besides 
some candidates of binar y MBHs (e.g..lSudou et al.ll2003l ; 
iBoroson fc LauCTi [200I 'P otti et all 120091 there is few 
evidence for multiple MBHs in a massive galaxy. Thus, 
multiple MBHs possibly merge into a heavier BH. But, 



the merger path of multiple MBHs in a massive galaxy 
has not been hitherto resolved. 

In this paper, we explore the evolution of the system 
of multiple MBHs in a massive galaxy. Previous stud- 
ies found that if two MBHs are in one galaxy, they are 
hard to merge within a Hubble time due to loss con e 
depletion (IBe gelman et al.lll980l; iMakino. Funatoll2004D . 
If three MBHs are in one galaxy, two of them merge 
or result in a binary, and the other is ejected from 
the galaxy (jlwasawa et al.l l2006l) . Here, we consider a 
two-component system that consists of MBHs and stars, 
where not only three-body interaction of MBHs but also 
dynamical friction by stars are incorporated. The pa- 
per is organized as follows. In section [21 we describe our 
simulation method to follow the evolution of the multiple 
MBHs. In section [3l simulation results are presented. In 
sectionlH we discuss the validity of our model. In section 
El we summarize this paper. 

2. METHOD 
2.1. Setup 

We treat MBHs and stars as a A^-body system. An 
individual MBH corresponds to one massive particle, and 
field stars composing a host galaxy are approximated as 
super particles, a fraction of which may be interpreted 
as dark matter. 

The field s tars are distributed according to the Hern- 
quist model ([Herngu ist 1990), whose radial mass density 
distribution, p(r), is given by 
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where M and are the total mass and virial radius of 
the galaxy, respectively. The virial radius is given by 
Tv = GM/4:\E\, where G is the gravitational constant, 
and E is the total energy of the galaxy. The number of 
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the field stars is iV = ^,12K{\K = 1024). 

Here, ten MBHs with equal mass of IO^Mq are set 
in a galaxy with 10^ ^M©. This means that the total 
mass of the MBHs is 0.1 % of the galaxy mass and the 
mass ratio of each MBH to each field star is about 50. 
We perform five simulations with different phase space 
distributions of the MBHs at the initial time in order to 
see the dependence on the stellar mass density around the 
MBHs. The MBHs are distributed initially within one- 
third, two-third, and one virial radius of the galaxy in 
model A, B, and C, respectively. Furthermore, in model 
A, the positions of MBHs are changed according to three 
different set of random number, which are labeled by 
model Ai, A2, and A3, respectively. These models are 
summarized in Table [T] 

In the present simulations, we adopt the standard A^- 
body units, G ~ M ^ ^ 1. In such units, the time 
unit of simulation, inu, is comparable to the dynamical 
time, tnu tdy The light speed is c = 600 in this units, 
which means that the three-dimensional velocity disper- 
sion of field stars is 300 km/s at the galaxy center. 

If we convert the A^-body units to physical units, the 
virial radius, the dynamical timescale within the virial 
radius, t^y, and the average mass density of the galaxy 
within the virial radius, p^, are respectively given by 
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According to this scaling, the present results can be ap- 
plicable for the different mass system, which is discussed 
later. 

2.2. Merger condition 

We assume that two MBHs merge through gravi- 
tational wave radiation, when the separation between 
two MBHs is less than ten times the sum of their 
Schwarzschild radii: 

|rB.i - rBjl < 10 (r-sch.i + '^schj) , (5) 

where tb.; and Tgch.i are the position and Schwarzschild 
radius of i-th MBH. The Schwarzschild radius of i-th 
MBH is Tsch.i = "^GmB.i/c^ for the MBH mass niB.i- 

2.3. Equation of motion 

The equations of motion for field stars and MBHs are 
respectively given by 
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where rf^.^ and rB,i are the positions of i-th field star 
and MBH, Ni and A'b are the numbers of field stars and 
MBHs, a.s,ij and afB,y are the accelerations of j-th field 



star and MBH on i-th field star, and SLBi.ij and aBB,ij 
are the accelerations of j-th field star and MBH on i-th 
MBH, respectively. Excepting the MBH-MBH interac- 
tion, the accelerations are given by Newtonian gravity: 
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where m{j and tob.j are respectively the masses of j-th 
field star and MBH, and the softening parameter (e = 
10"'^) is introduced only in star-star interactions. 

The acceleration between two MBHs contains Newto- 
nian gravity and post-Newtonian corrections, such as 
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For the second term (apN^i^), the pericenter 
(IPN and 2PN terms) as well as the gravita- 
tional i;adiB;ti£n__emissifln is consid- 
ered jDamour. fc DervelldHOSlI: ISoffelfTOSa IKupi et al.l 
20061). W e adop t equation (1), (2), (3), and (4) in 



Kupi et all (|2006 ) for the second term. 

If the semi-major axis is less than Ocrit = 5 x 10^^, 
the motion of the binary is transformed to the motion of 
the center of mass and the relative motion. We ignore 
tidal forces by distant field stars on the binary. Then, the 
acceleration by a distant field star k to the center of mass 
(acm.fc) and the relative motion (flrci,*;) is approximated 
as 



acm.A; ' 
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arcl,fc ~ 0, 



(12) 
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where rem is the position of the center of mass of the 
MBH binary. Distant field stars are defined as 



where C = 200. 



2.4. Numerical sche 
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We use a fourth-order Herrnite sc heme with individual 
timestep (iMakino fc AarsethI |1992|) and block timestep 
()McMillanlll986[) for field stars. MBHs, and the relative 
motion in binary MBHs. As for the motion of the center 
of mass of MB H binaries, a fourth-order H ermite Ahmad- 
Cohen scheme (jMakino fc Aarsetblll992l ) is employed. In 
an Hermite Ahmad- Cohen scheme, the acceleration due 
to single MBHs or stars near a binary MBH, and those 
due to distant stars are calculated on separate timesteps, 
which we call "neighbor step" and "distant step" , respec- 
tively. 

The timestep except the neighbor step is determined 

as 

At 



(15) 



where ry is the accuracy parameter, and a is the accelera- 
tion of a field star, a single MBH, or the relative motion 



Successive Merger of Massive Black Holes 



3 



TABLE 1 

Summary for our simulation results. 



Model name 


''MBh/^'v 


Rand # 




10*mB,scc 


Ai 


1/3 


Rl 


4 


3 (ejected) 


Aa 


1/3 


R2 


4 


1 


A3 


1/3 


R3 


6 


1 


B 


2/3 


Rl 


5 


1 


C 


1 


Rl 


5 


1 



of a binary MBH. The function / in equation ([15)) is 
given by 

|x||xW| + |x(^)p 
/W = i_.n^ii„rzt^i , I ^3)12 - (16) 



c(i)||x(4) 



The accuracy parameter is set to be 77 = 0.01 for timestep 
of field stars, 77 — 0.0025 for timestep of single MBHs and 
distant step of the center of mass of binary MBHs, and 
77 = 0.000625 for the relative motion in binary MBHs. 
The timestep for the center of mass of a binary MBH is 
determined as 

Ai = v/min[r;i/(a),772/(apN)], (17) 

where r/i = 0.0025 and 772 = 0.000625 are the accuracy 
parameters, and apN is the acceleration due to post- 
Newtonian corrections for the center of mass of a binary 
MBH. The relative error in our simulations is less than 
0.1 % regarding the total energy. 

We perform simulations with 64 nod es of the FIRST 
simul ator in University of Tsukuba (jUmemura et al.l 
120081) . At each node, the FIRST simulator is equipped 
with one Blade-GRAPE board, which is the accelerator 
of the gravity calculations for collisional A^-body prob- 
lem. The gravity by field stars for a given field star and 
MBH is calculated in parallel. 

3. RESULTS 
3.1. Merger of multiple MBHs 

We calculate a system of ten MBHs in one galaxy dur- 
ing about 140 A^-body time units, which corresponds to 
about 800 Myr in physical units. We find that several 
MBHs merge into one for all simulations. We summarize 
MBH masses at the final time of simulations in Table [TJ 
The heaviest MBH has TTie.max = 4 - 6 x 10^**, while 
the second heaviest MBH has 771b,scc = 1 x 10"**, ex- 
cept model Ai. The growth of such a dominant MBH 
is weakly dependent on the initial stellar mass density 
around MBHs. 

We see the process in which only one MBH grows in 
each simulation, using the result of model A3 as a typical 
case. As shown in the top panel of Figure [Jl only one 
MBH grows, and other MBHs do not grow. The second 
and third panels of Figure [T] show that MBHs compose a 
binary with a semi-major axis less than 10~* for a longer 
time than the dynamical time (idy '^1). When a binary 
merges to form a heavier MBH, a lighter component is 
often exchanged by a third MBH. Thereafter, a binary 
MBH forms again, containing the heaviest MBH (see the 
fourth panel of Figure [1]). This is because a heavier MBH 
is easier to be retained in a binary MBH through an 
interaction between the binary MBH and single MBH. 
Consequently, only one heavy MBH grows in a galaxy. 
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Fig. 1. — Time evolution of MBHs in model A3. From top to 
bottom, the panels show the mass of the heaviest MBH and second 
heaviest at each time, the semi-major axis, the distance at the 
pericenter, and the masses of two MBHs when a binary forms, and 
the number of binary MBHs. As for the number of binary MBHs, 
we only count binary MBHs whose semi-major axes are less than 
10~^. Pairs of integers in the second panel show the labels of 
MBHs composing the binary MBHs, where the heaviest MBH is 
labeled with " — 1". We attached labels only to binary MBHs which 
are long-lived, or merge eventually. In the second top and middle 
panels, filled circles indicate the moments when MBHs merge and 
crosses denote those when binary components are exchanged. 



The reason why no other MBHs merge is understood 
as follows. When a binary MBH forms in a galaxy, the 
binary prevents the formation of another binary MBH, 
because the preformed binary MBH gives other MBHs 
kick velocities. Since MBHs cannot merge without form- 
ing a binary MBH, no other MBHs can merge. 

In model Ai, the first binary forms just temporarily. It 
is ejected with a speed of more than 1000 km/s due to the 
back reaction of sling-shot mechanism, when the binary 
with masses of 2 x 10^** and 1 x 10^"* merges through the 
sling-shot mechanism of an interaction MBH. After the 
ejection, the second binary forms and merge to produce 
the MBH with mass of 4 x 10^"* through the same process 
as described above. In our simulations, the ejection by 
the sling-shot interaction is observed only once. In all 
the simulations, the merger of binary MBHs occurs 21 
times in total. Hence, the ejection of the heaviest MBH 
due to the sling-shot seems rare. 

3.2. Merger mechanism of a binary MBH 
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Fig. 2. — Time variation of the distance of the third nearest MBH 
from the galactic center. The result in model A3 and that in a fixed 
potential are compared. 

For the mergers of the MBHs, the dynamical friction 
plays a key role. In Figure [21 the time variation of the 
distance of the third nearest MBH from the galactic cen- 
ter is shown, where the result in the TV-body simulation 
and that in a fixed potential of the Hernquist model are 
compared. In the fixed potential, no dynamical friction 
is exerted. As seen in this figure, the dynamical friction 
by field stars in the iV-body simulation allows the MBHs 
to gather near the galaxy center. Thus, two MBHs can 
compose a binary MBH, and subsequently another MBH 
can intrude the binary MBH. In contrast, in a fixed po- 
tential, even a binary MBH can not be formed, since they 
can not gather at around the galaxy center. 

Here we see the details of the merger mechanism of 
binary MBHs. In Figure [3l the processes of the second 
merger of MBHs are shown in models A2 and A3. The 
second merger in model A2 is the simplest case. The bi- 
nary MBH and a single MBH approaches to each other, 
and strongly interact. Consequently, the distance of the 
binary MBH at the pericenter shrinks, followed by emit- 
ting gravitational wave radiation and merger. In prac- 
tice, another single MBH often intrudes into the binary 
MBH whose orbit is being decayed due to energy loss 
through gravitational wave radiation, and subsequently 
the semi-major axis and eccentricity of the binary MBH 
are changed, as seen in the middle and bottom panels of 
figure [3l However, the crucial impact is brought by one 
strong interaction. In our five simulations, there does 
not occur the simultaneous interaction of multiple MBHs 
that triggers the merger of a MBH binary. 

Also, we have not observed the secular angular mo- 
mentu m loss of a bi nary MBH through the Kozai mech- 
anism (|Kozaiiri962D . The Kozai mechanism can occur, 
only if the internal orbit of the binary MBH is closed 
in every binary period. However, the internal orbit is 
not closed du e to the relativistic pericenter shift (IPN 
and 2PN) (Bla es et al.ll2002l : iBerentzen et al.ll2009| ). Ac- 
tually, we have found that the Kozai mechanism works 
for a binary to merge only if we do not include the IPN 
and 2PN terms. The suppression of the Kozai mech- 
anism is als o demo nstrated in the case of stellar-sized 
black holes (iMiller fc Hamilton 2002) and in the plane- 
tary orbits (jFabrvckv fc Tremainell2007l) . 

4. DISCUSSION 

The present simulations have the scalings shown in 
equations ©, and Qj. Thus, the present results are 
applicable for the different mass scales of a host galaxy 
and MBHs, if the density of the host galaxy satisfies the 
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Fig. 3. — Time variation of the semi-major axes, a, and the 
distance at the pericenter of binary MBHs, rp, just before the 
second mergers in models A2 (top) and A3 (middle and bottom). 
Also, the separation between the center of mass of the binary MBHs 
and single MBHs which interact strongly with the binary MBHs is 
shown by dashed and dotted curves (rr, r4, ri, and r^), where the 
numbers attached to r are the labels of MBHs. The bottom panel 
demonstrates the exchange interaction that one binary component 
is replaced by another single MBH. 

scaling. For instance, if the collapse redshift of a host 
galaxy is shifted from 2; = 6toz=10, the virial density 
of the host galaxy is expected to be increased by a factor 
of 3.9. Hence, the same merger processes of MBHs are 
expected for a host galaxy with 2.6 x and MBHs 

with equal mass of 2.6 x lO^M©. 

In the present model, the host galaxy is assumed to 
have a spherical Hernquist profile. MBHs ejected from 
the center on a nearly radial orbit with speed below that 
of escape speed repeatedly return to the center in our 
spherical model. However, these MBHs do not contribute 
much to the merger of binary MBHs. We have found that 
there is no time for such MBHs to interact with binary 
MBHs, since they pass the center with high speed. MBHs 
which gather towards the galactic center by the dynam- 
ical friction more effectively contribute to the merger of 
a binary MBH. 

We have ignored the recoil by gravitational wave. 
Unless black hole spins are aligne d, the recoil veloc- 
ity can reach up t o 4000 km/s (Ca mpanelh et al.l 120071 : 
iLousto et al.ll201()h . Such a large recoil velocity can eject 
the merger remnant from the galaxy. However, if their 
spins are aligned b efore their merger due to relativis- 
tic spin precession ijKesden et al.|[20Tol) . then the recoil 
velocity decreases to a few 100 km/s. If the recoil ve- 
locity is a few lOOkm/s, the merger remnant is possibly 
confined in a host galaxy with the velocity dispersion of 
300 km/s. Then, the recoil may not the results dramat- 
ically. Nonetheless, it seems worth exploring carefully 
the effects of recoil, which will be considered in the fu- 
ture analysis. 

We have assumed that a binary MBH immediately 
merges at the moment when their separation become 
smaller than ten times the sum of their Schwarzschild 
radii. In practice, a binary MBHs take a bit more time 
to merge. But, the timescale is too short, the merger 
to be disturbed by the intrusion of another MBH. The 
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merging timescale of a binary MBH is ^ 10^^^ in TV-body 
units, while it interacts with a MBH once a dynamical 
time of the galaxy, 1 in iV-body units. Actually, as 
seen in Figure [U any MBH does not approach to the 
binary MBHs after the semi-major axis becomes about 
10"'' length unit. 

5. SUMMARY 

We have performed highly- accurate A^-body simula- 
tions to explore the evolution of multiple MBHs in one 
galaxy. Here, ten MBHs with equal mass of lO^Af© are 
set in a galaxy with 10^^ M©. As a result, it is found 
that 4-6 MBHs successively merge, resulting in a sin- 
gle heavier MBH within 1 Gyr. The growth timescale is 
shorter than a Hubble time at rcdshift z ^ 6. After 4 - 
6 MBHs merge, the other MBHs in the galaxy have the 
same masses as those at the initial time. 

The key physics for the successive merger is the dynam- 
ical friction that allows the formation of binary MBHs 
and frequent interactions of single MBHs with the bi- 
nary MBHs at the galactic center. Hence, the key pa- 
rameter for the MBH merger is the density of field stars 
in the regions where MBHs are distributed. The dis- 
tance of a MBH binary at the pericenter shrinks through 
the sling-shot mechanism of another MBH. It followed 
by the shrink of the semi-major axis due to the energy 
loss by gravitational wave radiation, and eventually the 
binary merges. We have found that the secular angular 
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is properly included. The present simulations imply that 
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cess. The present results are applicable for the different 
mass scales of a host galaxy and MBHs, if the density of 
the host galaxy satisfies the scaling in the A^-body units. 
This merger path can be important for the growth of 
MBHs in a primordial galaxy. Also, the MBH merger 
may be a potential source of gravitational waves for the 
Laser Interferometer Space Antenna (LISA) and pulsar 
timing. 
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